#include <cstdio>
#include <cstdlib>
#include <string>
#include <omp.h>
#include <cmath>
#ifdef _WIN64
#include "windows.h"
#else
#include <time.h>
#endif
#include "limits.h"
#include "float.h"
#include <vector>

static double* a = nullptr;
static double* b = nullptr;
static double* c = nullptr;
static int N = 0;
static double	avgtime[4] = { 0 }, maxtime[4] = { 0 }, mintime[4] = { FLT_MAX,FLT_MAX,FLT_MAX,FLT_MAX };
static char* label[4] = { "Copy:      ", "Scale:     ", "Add:       ", "Triad : " };
# define NTIMES	10
# define OFFSET	0

# ifndef MIN
# define MIN(x,y) ((x)<(y)?(x):(y))
# endif
# ifndef MAX
# define MAX(x,y) ((x)>(y)?(x):(y))
# endif

# define HLINE "-------------------------------------------------------------\n"

static void checkSTREAMresults()
{
	double aj, bj, cj, scalar;
	double asum, bsum, csum;
	double epsilon;
	int	j, k;

	/* reproduce initialization */
	aj = 1.0;
	bj = 2.0;
	cj = 0.0;
	/* a[] is modified during timing check */
	aj = 2.0E0 * aj;
	/* now execute timing loop */
	scalar = 3.0;
	for (k = 0; k < NTIMES; k++)
	{
		cj = aj;
		bj = scalar * cj;
		cj = aj + bj;
		aj = bj + scalar * cj;
	}
	aj = aj * (double)(N);
	bj = bj * (double)(N);
	cj = cj * (double)(N);

	asum = 0.0;
	bsum = 0.0;
	csum = 0.0;
	for (j = 0; j < N; j++) {
		asum += a[j];
		bsum += b[j];
		csum += c[j];
	}
#ifdef VERBOSE
	printf("Results Comparison: \n");
	printf("        Expected  : %f %f %f \n", aj, bj, cj);
	printf("        Observed  : %f %f %f \n", asum, bsum, csum);
#endif

#define abs(a) ((a) >= 0 ? (a) : -(a))
	epsilon = 1.e-8;

	if (abs(aj - asum) / asum > epsilon) {
		printf("Failed Validation on array a[]\n");
		printf("        Expected  : %f \n", aj);
		printf("        Observed  : %f \n", asum);
	}
	else if (abs(bj - bsum) / bsum > epsilon) {
		printf("Failed Validation on array b[]\n");
		printf("        Expected  : %f \n", bj);
		printf("        Observed  : %f \n", bsum);
	}
	else if (abs(cj - csum) / csum > epsilon) {
		printf("Failed Validation on array c[]\n");
		printf("        Expected  : %f \n", cj);
		printf("        Observed  : %f \n", csum);
	}
	else {
		printf("Solution Validates\n");
	}
}

double mysecond() {
#ifdef _WIN64
	//用QueryPerformanceCounter()来计时  微秒
	LARGE_INTEGER  large_interger;
	double dff;
	__int64  c1;
	QueryPerformanceFrequency(&large_interger);
	dff = large_interger.QuadPart;
	QueryPerformanceCounter(&large_interger);
	c1 = large_interger.QuadPart;
	return (double)(c1) / dff;
#else
	struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return 1. * t.tv_sec + 1.e-9 * t.tv_nsec;
#endif
}

#ifdef _WIN64
void test_mysecond()
{
	//用QueryPerformanceCounter()来计时  微秒
	LARGE_INTEGER  large_interger;
	double dff;
	__int64  c1, c2;
	QueryPerformanceFrequency(&large_interger);
	dff = large_interger.QuadPart;
	QueryPerformanceCounter(&large_interger);
	c1 = large_interger.QuadPart;
	Sleep(800);
	QueryPerformanceCounter(&large_interger);
	c2 = large_interger.QuadPart;
	printf("Counter freq on this machine %lf\n", dff);
	printf("First counter %I64d second counter %I64d diff %I64d\n", c1, c2, c2 - c1);
	printf("Elapsed % lf ms\n", (c2 - c1) * 1000 / dff);
}
#endif

# define M 20

int checktick()
{
	int		i, minDelta, Delta;
	double	t1, t2, timesfound[M];

	/*  Collect a sequence of M unique time values from the system. */

	for (i = 0; i < M; i++) {
		t1 = mysecond();
		while (((t2 = mysecond()) - t1) < 1.0E-6)
			;
		timesfound[i] = t1 = t2;
	}

	/*
	 * Determine the minimum difference between these M values.
	 * This result will be our estimate (in microseconds) for the
	 * clock granularity.
	 */

	minDelta = 1000000;
	for (i = 1; i < M; i++) {
		Delta = (int)(1.0E6 * (timesfound[i] - timesfound[i - 1]));
		minDelta = MIN(minDelta, MAX(Delta, 0));
	}

	return(minDelta);
}

void stream_test(int Nelem)
{
#ifdef _WIN64
	test_mysecond();
#endif
	int			quantum;
	int			BytesPerWord;
	int	j, k;
	double		scalar, t, times[4][NTIMES];
	printf("Number of Threads requested = %i\n", omp_get_max_threads());
	/* --- SETUP --- determine precision and check timing --- */

	printf(HLINE);
	BytesPerWord = sizeof(double);
	printf("This system uses %d bytes per DOUBLE PRECISION word.\n",
		BytesPerWord);

	N = Nelem;
	a = new double[N];
	b = new double[N];
	c = new double[N];
	const double bytes[4] = {
		(double)2 * sizeof(double) * N,
		(double)2 * sizeof(double) * N,
		(double)3 * sizeof(double) * N,
		(double)3 * sizeof(double) * N
	};

	printf(HLINE);
	printf("Array size = %d, Offset = %d\n", N, OFFSET);
	printf("Total memory required = %.1f MB.\n",
		(3.0 * BytesPerWord) * ((double)N / 1048576.0));
	printf("Each test is run %d times, but only\n", NTIMES);
	printf("the *best* time for each is used.\n");
#ifdef _WIN64
	auto getMask =
		[](HANDLE handle)
	{
		auto mask = DWORD_PTR{ 1 };
		mask = SetThreadAffinityMask(handle, mask);
		SetThreadAffinityMask(handle, mask);
		return mask;
	};
	std::vector< DWORD_PTR > masks(omp_get_max_threads());
	#pragma omp parallel
	{
		int tid = omp_get_thread_num();
		masks[tid] = getMask(GetCurrentThread());
		printf("tid %d Mask %d\n", tid, masks[tid]);
	}
#endif
	/* Get initial value for system clock. */
	#pragma omp parallel for schedule(static)
	for (j = 0; j < N; j++) {
		a[j] = 1.0;
		b[j] = 2.0;
		c[j] = 0.0;
	}

	printf(HLINE);

	if ((quantum = checktick()) >= 1)
		printf("Your clock granularity/precision appears to be "
			"%d microseconds.\n", quantum);
	else
		printf("Your clock granularity appears to be "
			"less than one microsecond.\n");

	t = mysecond();

	if ((quantum = checktick()) >= 1)
		printf("Your clock granularity/precision appears to be %d microseconds.\n", quantum);
	else
		printf("Your clock granularity appears to be less than one microsecond.\n");

	t = mysecond();
	
	#pragma omp parallel for schedule(static)
	for (j = 0; j < N; j++)
		a[j] = 2.0E0 * a[j];
	t = 1.0E6 * (mysecond() - t);

	printf("Each test below will take on the order of %d microseconds.\n", (int)t);
	printf("   (= %d clock ticks)\n", (int)(t / quantum));
	printf("Increase the size of the arrays if this shows that\n");
	printf("you are not getting at least 20 clock ticks per test.\n");

	printf(HLINE);

	printf("WARNING -- The above is only a rough guideline.\n");
	printf("For best results, please be sure you know the\n");
	printf("precision of your system timer.\n");
	printf(HLINE);

	/*	--- MAIN LOOP --- repeat test cases NTIMES times --- */
	scalar = 3.0;
	for (k = 0; k < NTIMES; k++)
	{
		times[0][k] = mysecond();

		#pragma omp parallel for schedule(static)
		for (j = 0; j < N; j++)
			c[j] = a[j];

		times[0][k] = mysecond() - times[0][k];

		times[1][k] = mysecond();

		#pragma omp parallel for schedule(static)
		for (j = 0; j < N; j++)
			b[j] = scalar * c[j];

		times[1][k] = mysecond() - times[1][k];

		times[2][k] = mysecond();

		#pragma omp parallel for schedule(static)
		for (j = 0; j < N; j++)
			c[j] = a[j] + b[j];

		times[2][k] = mysecond() - times[2][k];

		times[3][k] = mysecond();

		#pragma omp parallel for schedule(static)
		for (j = 0; j < N; j++)
			a[j] = b[j] + scalar * c[j];

		times[3][k] = mysecond() - times[3][k];
	}

	/*	--- SUMMARY --- */

	for (k = 1; k < NTIMES; k++) {/* note -- skip first iteration */
		for (j = 0; j < 4; j++) {
			avgtime[j] = avgtime[j] + times[j][k];
			mintime[j] = MIN(mintime[j], times[j][k]);
			maxtime[j] = MAX(maxtime[j], times[j][k]);
		}
	}

	printf("Function      Rate (MB/s)   Avg time     Min time     Max time\n");
	for (j = 0; j < 4; j++) {
		avgtime[j] = avgtime[j] / (double)(NTIMES - 1);

		printf("%s%11.4f  %11.4f  %11.4f  %11.4f\n", label[j],
			1.0E-06 * bytes[j] / mintime[j],
			avgtime[j],
			mintime[j],
			maxtime[j]);
	}
	printf(HLINE);

	/* --- Check Results --- */
	checkSTREAMresults();
	printf(HLINE);

	delete[] a; a = nullptr;
	delete[] b; b = nullptr;
	delete[] c; c = nullptr;
}

int main(int argc, char* argv[])
{
    int N = atoi(argv[1]);
    stream_test(N);

    return 0;
}